\section{方差缩减方法}

随机模拟方法虽然有着适用性广、方法简单的优点, 但是又有精度低、计算量大的缺点, 一整套模拟算几天几夜也是常有的事情。如果能成倍地减小随机模拟误差方差, 就可以有效地节省随机模拟时间, 有些情况下可以把耗时长到不具有可行性的模拟计算 (比如几个月) 缩短到可行 (比如几天)。

{\color{red}\textbf{[为什么方差缩减如此重要]:} 随机模拟的误差与$1/\sqrt{N}$成正比。如果方差缩减10倍,就相当于样本量增加100倍——这可能意味着从几个月缩短到几天。方差缩减是``用智慧替代算力''的典型例子。}

前一节的关于定积分计算的重要抽样法、分层抽样法都是降低随机模拟误差方差的重要方法, 也可以用在一般的模拟问题中。本节介绍一些其它的方差缩减技巧。我们以随机变量 \(X\) 的期望 \(\theta  = {EX}\) 的估计为例,目标是降低 \(\theta\) 的估计量的渐近方差。

{\color{red}\textbf{[方差缩减的四大武器]:} 控制变量法——利用已知信息校正偏差;对立变量法——利用对称性产生负相关;条件期望法——用解析公式消除部分随机性;随机数复用——配对比较时误差抵消。这些方法可单独使用,有些可以组合,但核心思想都是``减少独立随机源的数量''。}

\subsection{控制变量法}

设要估计随机变量 \(X\) 的期望 \(\theta  = {EX}\) ,从 \(X\) 中抽取 \(N\) 个独立样本值 \({X}_{1},{X}_{2},\ldots ,{X}_{n}\) , 用样本平均值 \(\bar{X} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{X}_{i}\) 估计 \({EX}\) 。为了提高精度可以利用辅助信息。设有另外的随机变量 \(Y\) 满足

\[
{EY} = 0,\;\operatorname{Cov}\left( {X,Y}\right)  < 0
\]

{\color{red}\textbf{[核心思想]:} 控制变量法的本质是``借用已知信息''。如果$Y$与$X$负相关且期望为零,那么当$X$高估时$Y$倾向于为正(可以减去),当$X$低估时$Y$倾向为负(相当于加上)。这种``自动校正''机制能显著降低方差。}

令 \(Z = X + Y\) ,则

\[
E\left( Z\right)  = \theta ,\;\operatorname{Var}\left( Z\right)  = \operatorname{Var}\left( X\right)  + \operatorname{Var}\left( Y\right)  + 2\operatorname{Cov}\left( {X,Y}\right) ,
\]

只要 \(\operatorname{Var}\left( Y\right)  + 2\operatorname{Cov}\left( {X,Y}\right)  < 0\) 则 \(\operatorname{Var}\left( Z\right)  < \operatorname{Var}\left( X\right)\) ,如果有(X, Y)成对的抽样 \(\left( {{X}_{i},{Y}_{i}}\right) ,i =\)  \(1,2,\ldots ,n\) ,令 \({Z}_{i} = {X}_{i} + {Y}_{i}\) ,则用 \(\bar{Z}\) 来估计 \(\theta  = {EX} = {EZ}\) 的渐近方差就比用 \(\bar{X}\) 估计 \(I\) 的渐近方差减小了。

为了最好地利用 \(Y\) 与 \(X\) 的相关性,令

\[
Z\left( b\right)  = X + {bY},
\]

则

\[
{EZ}\left( b\right)  = {EX} = \theta ,
\]

\[
\operatorname{Var}\left( {Z\left( b\right) }\right)  = \operatorname{Var}\left( X\right)  + {2b}\operatorname{Cov}\left( {X,Y}\right)  + {b}^{2}\operatorname{Var}\left( Y\right) ,
\]

求 \(\operatorname{Var}\left( {Z\left( b\right) }\right)\) 关于 \(b\) 的最小值点,得

\[
b =  - \operatorname{Cov}\left( {X,Y}\right) /\operatorname{Var}\left( Y\right)  =  - {\rho }_{X,Y}\sqrt{\operatorname{Var}\left( X\right) /\operatorname{Var}\left( Y\right) },
\]

这时

\[
\operatorname{Var}\left( {Z\left( b\right) }\right)  = \left( {1 - {\rho }_{X,Y}^{2}}\right) \operatorname{Var}\left( X\right)  \leq  \operatorname{Var}\left( X\right) ,
\]

{\color{red}\textbf{[关键公式]:} 这个$(1-\rho_{X,Y}^2)$因子揭示了方差缩减的威力——如果相关系数$|\rho_{X,Y}|=0.99$(强相关),方差就降至原来的$1-0.99^2=0.02$,即缩减50倍!这就是为什么寻找高度相关的控制变量如此重要。}

可见只要能找到零均值随机变量 \(Y\) 使得 \({\rho }_{X,Y} \neq  0\) 就可以减小 \({EX}\) 的估计方差。 \(Y\) 和 \(X\) 的相关性越强, 改善幅度越大。这种减小随机模拟误差方差的方法叫做控制变量法。

实际中 \({\rho }_{X,Y}\) 和 \(\operatorname{Var}\left( X\right) ,\operatorname{Var}\left( Y\right)\) 可能是未知的,可以先模拟一个小的样本估计 \({\rho }_{X,Y}\) 和 \(\operatorname{Var}\left( X\right) ,\operatorname{Var}\left( Y\right)\) 从而获得 \(b\) 的估计值。

控制变量法中要求控制变量 \(Y\) 与 \(X\) 相关且 \({EY} = 0\) 。如果 \({EY} \neq  0\) 但 \({EY} = {\mu }_{Y}\) 已知, 只要用 \(Y - {\mu }_{Y}\) 代替 \(Y\) ,这需要能预先知道 \(Y\) 的期望值的真值。另一种情况是 \({EY} = {EX}\) 未知, \(Y\) 与 \(X\) 相关,这时令 \(Z = {\alpha X} + \left( {1 - \alpha }\right) Y,\alpha  \in  \left\lbrack  {0,1}\right\rbrack\) ,则 \({EZ} = {EX} = \theta\) ,可以求 \(\alpha\) 使得 \(\operatorname{Var}\left( Z\right)\) 最小。容易知道当 \(\alpha  = \operatorname{Cov}\left( {Y,Y - X}\right) /\operatorname{Var}\left( {Y - X}\right)\) 时 \(\operatorname{Var}\left( Z\right)\) 最小。

例 3.3.1. 设要估计 \(I = {\int }_{0}^{1}{e}^{t}{dt}\) 。当然,可以得到积分真值为 \(e - 1\) ,这里用来演示控制变量法的优势。

{\color{red}\textbf{[例子策略]:} 虽然真值已知,但这个简单例子能清晰展示方差缩减效果——可以精确计算理论方差并与模拟对比。复杂问题中无法算出理论方差,但缩减原理相同。}

设 \(U \sim  \mathrm{U}\left( {0,1}\right) ,X = {e}^{U}\) ,则 \(I = E{e}^{U} = {EX}\) ,可以用平均值法估计 \(\mathrm{I}\) 为

\[
{\widehat{I}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{e}^{{U}_{i}}. \tag{3.45}
\]

其方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  = \frac{1}{N}\operatorname{Var}\left( {e}^{U}\right)  = \frac{1}{N}\left( {-\frac{1}{2}{e}^{2} + {2e} - \frac{3}{2}}\right)  \approx  \frac{0.2420}{N}. \tag{3.46}
\]

令 \(Y = U - \frac{1}{2}\) ,则 \({EY} = 0,X\) 与 \(Y\) 正相关,可以计算出 \(\operatorname{Cov}\left( {X,Y}\right)  \approx  {0.14086},\operatorname{Var}\left( Y\right)  = 1/{12}\) (更复杂的问题中可能需要从一个小的随机抽样中近似估计),于是 \(b =  - \operatorname{Cov}\left( {X,Y}\right) /\operatorname{Var}\left( Y\right)  =\) -1.690,对 \(Z\left( b\right)  = {e}^{U} - {1.690}\left( {U - \frac{1}{2}}\right)\) 有 \(\operatorname{Var}\left( {Z\left( b\right) }\right)  = \left\lbrack  {1 - {\rho }_{X,Y}^{2}}\right\rbrack  \operatorname{Var}\left( X\right)  = \left( {1 - {0.9919}^{2}}\right) \operatorname{Var}\left( X\right)  =\)  \({0.016}\operatorname{Var}\left( X\right)  = {0.0039}\) ,用控制变量法估计 \(I\) 为

\[
{\widehat{I}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\left\lbrack  {{e}^{{U}_{i}} - {1.690}\left( {{U}_{i} - \frac{1}{2}}\right) }\right\rbrack  .
\]

\({\widehat{I}}_{1}\) 的方差比控制变量法 \({\widehat{I}}_{2}\) 的方差大 60 倍以上。

{\color{red}\textbf{[量级说明]:} 60倍方差缩减意味着达到相同精度只需$1/60$的样本量,或者说相同样本下误差标准差降至原来的$1/\sqrt{60}\approx 1/8$。这就是为什么控制变量法在金融蒙特卡罗中如此重要——能把几小时的计算缩短到几分钟。}

例 3.3.2. (系统可靠性估计) 考虑由 \(n\) 个部件组成的一个系统,用 \({S}_{i}\) 表示第 \(i\) 个部件是否正常工作,1 表示正常工作,0 表示失效。设 \({S}_{i} \sim  \mathrm{B}\left( {1,{p}_{i}}\right)\) 且各 \({S}_{i}\) 相互独立。用 \(Y\) 表示系统是否工作正常,1 表示工作正常,0 表示系统失效。设 \(Y\) 为 \({S}_{1},{S}_{2},\ldots ,{S}_{n}\) 的函数 \(\phi \left( {{S}_{1},{S}_{2},\ldots ,{S}_{n}}\right)\) 且 \(\phi\) 关于每个 \({S}_{i}\) 是单调不减,称 \(\phi\) 为系统的结构函数。令 \(R = P\left( {Y = 1}\right)  = {E\phi }\left( {{S}_{1},{S}_{2},\ldots ,{S}_{n}}\right)\) , 称 \(R\) 为系统可靠度。

例如, \(\phi \left( {{s}_{1},{s}_{2},\ldots ,{s}_{n}}\right)  = \mathop{\prod }\limits_{{i = 1}}^{n}{s}_{i}\) ,则当且仅当所有部件正常工作时系统才正常工作,这样的系统称为串联系统, 这时系统可靠度为

\[
R = P\left( {{S}_{1} = 1,{S}_{2} = 1,\ldots ,{S}_{n} = 1}\right)  = \mathop{\prod }\limits_{{i = 1}}^{n}P\left( {{S}_{i} = 1}\right)  = {p}_{1}{p}_{2}\ldots {p}_{n}. \tag{3.47}
\]

在系统比较简单的情况下 (如串联、并联),可以给出用 \({p}_{1},{p}_{2},\ldots ,{p}_{n}\) 表示 \(R\) 的表达式。但是更复杂的系统则很难写出 \(R\) 的表达式,这时可以用随机模拟方法估计 \(R\) 。记 \(\mathbf{S} = \left( {{S}_{1},{S}_{2},\ldots ,{S}_{n}}\right) ,\mathbf{X} = \phi \left( \mathbf{S}\right)\) ,对 \(\mathbf{S}\) 独立抽取 \(N\) 个点 \({\mathbf{S}}^{\left( j\right) } = \left( {{S}_{1}^{\left( j\right) },{S}_{2}^{\left( j\right) },\ldots ,{S}_{n}^{\left( j\right) }}\right) ,j =\)  \(1,2,\ldots ,N,R\) 可以用平均值法估计为

\[
{\widehat{R}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}\phi \left( {\mathbf{S}}^{\left( j\right) }\right) . \tag{3.48}
\]

令 \(Y = \mathop{\sum }\limits_{{i = 1}}^{n}\left( {{S}_{i} - {p}_{i}}\right)\) ,则 \({EY} = 0,Y\) 与 \(X\) 正相关。用一个小的抽样先近似估计 \(\operatorname{Cov}\left( {X,Y}\right) ,\operatorname{Var}\left( Y\right)\) 得到 \(b\) 的近似值,可以得到方差缩减的估计量

\[
{\widehat{R}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}\left\lbrack  {\phi \left( {\mathbf{S}}^{\left( j\right) }\right)  + b\mathop{\sum }\limits_{{i = 1}}^{n}\left( {{S}_{i}^{\left( j\right) } - {p}_{i}}\right) }\right\rbrack  .
\]

\subsection{对立变量法}

控制变量法需要知道控制变量 \(Y\) 的期望值真值,并精确或近似知道 \(\operatorname{Cov}\left( {X,Y}\right)\) 和 \(\operatorname{Var}\left( Y\right)\) 。对立变量法的要求比较简单。

{\color{red}\textbf{[对比优势]:} 控制变量法需要计算协方差和方差——有时很困难。对立变量法更巧妙:利用对称性自动产生负相关,无需计算任何统计量!这是``几何对称性''战胜``代数计算''的典型例子。}

模拟中经常使用均匀分布随机数 \(U\) 变换产生的随机数 \(X = g\left( U\right)\) 。下面的定理说明,如果变换 \(g\left( \cdot \right)\) 是单调的,则随机变量 \(Y = g\left( {1 - U}\right)\) 就是与 \(X\) 负相关的。注意 \(g\left( {1 - U}\right)\) 与 \(g\left( U\right)\) 同分布所以 \({EY} = {EX}\) 。

{\color{red}\textbf{[几何直觉]:} 想象$[0,1]$区间上的``镜像对称''——如果$U$在左边($U<0.5$),则$1-U$在右边;如果$U$偏大,则$1-U$偏小。对单调函数$g$,一个高另一个就低,天然负相关!}

定理 3.3.1. 设 \(g\) 为单调函数, \(U \sim  \mathrm{U}\left( {0,1}\right)\) ,则 \(\operatorname{Cov}\left( {g\left( U\right) ,g\left( {1 - U}\right) }\right)  \leq  0\) 。

证明. \(\forall {u}_{1},{u}_{2} \in  \left\lbrack  {0,1}\right\rbrack\) ,由 \(g\) 单调可知

\[
\left( {g\left( {u}_{1}\right)  - g\left( {u}_{2}\right) }\right) \left( {g\left( {1 - {u}_{1}}\right)  - g\left( {1 - {u}_{2}}\right) }\right)  \leq  0
\]

设 \({U}_{2}\) 服从 \(\mathrm{U}\left( {0,1}\right)\) 且与 \(U\) 独立,令 \({X}_{1} = g\left( U\right) ,{Y}_{1} = g\left( {1 - U}\right) ,{X}_{2} = g\left( {U}_{2}\right) ,{Y}_{2} = g\left( {1 - {U}_{2}}\right)\) , 则 \({X}_{1},{Y}_{1},{X}_{2},{Y}_{2}\) 的分布相同,且

\[
E\left( {{X}_{1} - {X}_{2}}\right) \left( {{Y}_{1} - {Y}_{2}}\right)
\]

\[
= \operatorname{Cov}\left( {{X}_{1} - {X}_{2},{Y}_{1} - {Y}_{2}}\right)
\]

\[
= \operatorname{Cov}\left( {{X}_{1},{Y}_{1}}\right)  + \operatorname{Cov}\left( {{X}_{2},{Y}_{2}}\right)  - \operatorname{Cov}\left( {{X}_{1},{Y}_{2}}\right)  - \operatorname{Cov}\left( {{X}_{2},{Y}_{1}}\right)
\]

\[
= 2\operatorname{Cov}\left( {{X}_{1},{Y}_{1}}\right)
\]

注意 \(\left( {{X}_{1} - {X}_{2}}\right) \left( {{Y}_{1} - {Y}_{2}}\right)  \leq  0\) 所以 \(\operatorname{Cov}\left( {{X}_{1},{Y}_{1}}\right)  \leq  0\) ,即 \(\operatorname{Cov}\left( {g\left( U\right) ,g\left( {1 - U}\right) }\right)  \leq  0\) 。证毕。 定理3.3.1可以推广到如下情形。

定理 3.3.2. 设 \(h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 是关于每个自变量单调的函数, \({U}_{1},{U}_{2},\ldots ,{U}_{n}\) 相互独立,则 \(\operatorname{Cov}\left( {h\left( {{U}_{1},{U}_{2},\ldots ,{U}_{n}}\right) ,h\left( {1 - {U}_{1},1 - {U}_{2},\ldots ,1 - {U}_{n}}\right) }\right)  \leq  0\) 。

证明略去,参见 Ross \({\left( {2013}\right) }^{\left\lbrack  {34}\right\rbrack  }§{9.9}\) 。

对均匀随机数 \(U\) 最常见的变换是逆变换 \(X = {F}^{-1}\left( U\right)\) 。下面的定理给出了提高 \(I = {EX}\) 估计精度的方法。

定理 3.3.3 (对立变量法). 设 \(F\left( x\right)\) 为连续分布函数, \(U \sim  \mathrm{U}\left( {0,1}\right) ,X = {F}^{-1}\left( U\right) ,Y =\)  \({F}^{-1}\left( {1 - U}\right) ,Z = \frac{X + Y}{2}\) ,则 \(X\) 与 \(Y\) 同分布 \(F\left( x\right)\) 且 \(\operatorname{Cov}\left( {X,Y}\right)  \leq  0\) ,

\[
\operatorname{Var}\left( Z\right)  \leq  \frac{1}{2}\operatorname{Var}\left( X\right)
\]

{\color{red}\textbf{[白捡的好处]:} 这个定理说:只需对每个$U_i$同时计算$F^{-1}(U_i)$和$F^{-1}(1-U_i)$取平均,方差就至少减半!计算量几乎不变(只多一次函数求值),效果却等价于样本量翻倍。这是最``划算''的方差缩减技巧。}

证明. 因为 \(U\) 和 \(1 - U\) 同分布所以 \(X = {F}^{-1}\left( U\right)\) 和 \(Y = {F}^{-1}\left( {1 - U}\right)\) 同分布。由定理3.3.1, 令 \(g\left( \cdot \right)  = {F}^{-1}\left( \cdot \right)\) 可知 \(\operatorname{Cov}\left( {X,Y}\right)  = \operatorname{Cov}\left( {g\left( U\right) ,g\left( {1 - U}\right) }\right)  \leq  0\) ,从而

\[
\operatorname{Var}\left( Z\right)  = \frac{\operatorname{Var}\left( X\right)  + \operatorname{Var}\left( Y\right)  + 2\operatorname{Cov}\left( {X,Y}\right) }{4}
\]

\[
= \frac{\operatorname{Var}\left( X\right)  + \operatorname{Cov}\left( {X,Y}\right) }{2} \leq  \frac{1}{2}\operatorname{Var}\left( X\right)
\]

证毕。

根据定理 3.3.3 的结论,为了估计 \(I = {EX}\) ,产生 \({U}_{1},{U}_{2},\ldots ,{U}_{N}\) 后用

\[
{Z}_{i} = \frac{1}{2}\left( {{F}^{-1}\left( {U}_{i}\right)  + {F}^{-1}\left( {1 - {U}_{i}}\right) }\right) ,i = 1,2,\ldots ,N
\]

的样本平均值估计 \(I\) 可以提高精度,在不增加抽样个数的条件下把估计的随机误差方差降低到原来的 \(\frac{1}{2}\) 以下。这种提高随机模拟精度的方法叫做对立变量法。对立变量法不需要计算 \(X\) 和 \(Y\) 的方差及协方差的值,比控制变量法更简便易行。

例 3.3.3. 再次考虑例3.3.1的问题,估计 \(I = {\int }_{0}^{1}{e}^{t}{dt}\) 。下面用对立变量法改善原始的平均值法 \({\widehat{I}}_{1}\) 的估计方差。

设 \(U \sim  \mathrm{U}\left( {0,1}\right) ,X = {e}^{U}\) ,令 \(Y = {e}^{1 - U}\) ,用对立变量法估计 \(I\) 为

\[
{\widehat{I}}_{3} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{{e}^{{U}_{i}} + {e}^{1 - {U}_{i}}}{2}, \tag{3.49}
\]

方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{3}\right)  = \frac{1}{N}\frac{\operatorname{Var}\left( {e}^{U}\right)  + \operatorname{Cov}\left( {{e}^{U},{e}^{1 - U}}\right) }{2} = \frac{1}{N}\left( {-\frac{3}{4}{e}^{2} + \frac{5}{2}e - \frac{5}{4}}\right)  \approx  \frac{0.003913}{N}, \tag{3.50}
\]

\({\widehat{I}}_{1}\) 的方差比对立变量法估计 \({\widehat{I}}_{3}\) 的方差大至少 60 倍,而 \({\widehat{I}}_{3}\) 的方差和例3.3.1中控制变量法估计量 \({\widehat{I}}_{2}\) 的方差相近。

{\color{red}\textbf{[惊人巧合]:} 对立变量法不需要计算协方差,却达到了与精心设计的控制变量法相近的效果!这是因为$U$和$1-U$的对称性恰好产生了接近最优的负相关。这种``免费午餐''在许多问题中都存在。}

对立变量法和控制变量法是类似做法, 一般不能结合使用。

{\color{red}\textbf{[使用场景]:} 控制变量法适合有``参考量''的情况(如金融中的Delta对冲);对立变量法适合逆变换法生成随机数的情况;条件期望法适合分层模型;随机数复用适合算法比较。实践中先尝试最简单的对立变量法,效果不佳再考虑控制变量法。}

例 3.3.4. 再次考虑例3.3.2的可靠度估计问题。用对立变量法改善估计方差。设 \(\left\{  {U}_{k}\right\}\) 为标准均匀分布随机数列, 取

\[
{S}_{i}^{\left( j\right) } = \left\{  \begin{array}{ll} 1 & \text{ 当 }{U}_{n\left( {j - 1}\right)  + i} \leq  {p}_{i}, \\  0 & \text{ 其它 } \end{array}\right.  \tag{3.51}
\]

则 \({S}_{i}^{\left( j\right) }\) 是 \({U}_{n\left( {j - 1}\right)  + i}\) 的单调非增函数。 \(R\) 用平均值法估计为

\[
{\widehat{R}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}\phi \left( {{S}_{1}^{\left( j\right) },{S}_{2}^{\left( j\right) },\ldots ,{S}_{n}^{\left( j\right) }}\right) . \tag{3.52}
\]

利用对立变量法,令 \(h\left( {{U}_{1},{U}_{2},\ldots ,{U}_{n}}\right)  = \phi \left( {{S}_{1},{S}_{2},\ldots ,{S}_{n}}\right)\) ,则 \(h\) 关于每个自变量是单调非增函数, 于是

\[
\operatorname{Cov}\left( {h\left( {{U}_{1},{U}_{2},\ldots ,{U}_{n}}\right) ,h\left( {1 - {U}_{1},1 - {U}_{2},\ldots ,1 - {U}_{n}}\right) }\right)  \leq  0, \tag{3.53}
\]

估计系统可靠度 \(R\) 为

\[
{\widehat{R}}_{3} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}\frac{h\left( {{U}_{1}^{\left( j\right) },{U}_{2}^{\left( j\right) },\ldots ,{U}_{n}^{\left( j\right) }}\right)  + h\left( {1 - {U}_{1}^{\left( j\right) },1 - {U}_{2}^{\left( j\right) },\ldots ,1 - {U}_{n}^{\left( j\right) }}\right) }{2} \tag{3.54}
\]

就能比 \({\widehat{R}}_{1}\) 的误差方差至少降低 \(\frac{1}{2}\) 。

注意定理3.3.2中的 \({U}_{1},{U}_{2},\ldots ,{U}_{n}\) 仅要求独立,并未指定分布。事实上,此结果还可以推广。如果 \(\left( {{X}_{i},{Y}_{i}}\right) ,i = 1,2,\ldots ,n\) 独立,且 \({X}_{i}\) 和 \({Y}_{i}\) 负相关, \(h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 关于每个自变量单调不减,则 \(\operatorname{Cov}\left( {h\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right) ,h\left( {{Y}_{1},{Y}_{2},\ldots ,{Y}_{n}}\right) }\right)  \leq  0\) ,可以类似地用控制变量法或对立变量法的做法构造方差缩减估计量。例如,设 \({X}_{i} \sim  \mathrm{N}\left( {{\mu }_{i},{\sigma }_{i}^{2}}\right)\) 相互独立,要估计 \(\theta  =\)  \({Eh}\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right)\) ,其中 \(h\) 关于每个自变量单调不减,则 \(h\left( {2{\mu }_{1} - {X}_{1},2{\mu }_{2} - {X}_{2},\ldots ,2{\mu }_{n} - {X}_{n}}\right)\) 与 \(h\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right)\) 同分布,对

\[
Z = \frac{h\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right)  + h\left( {2{\mu }_{1} - {X}_{1},2{\mu }_{2} - {X}_{2},\ldots ,2{\mu }_{n} - {X}_{n}}\right) }{2}
\]

抽样用平均值估计 \(\theta\) 可以比仅对 \(h\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right)\) 抽样得到的估计量的方差缩减一半以上。

\subsection{条件期望法}

进行统计估计时, 如果有额外的相关信息, 利用这样的信息可以提高估计精度。比如, 对随机变量 \(Y\) ,如果 \(Y\) 服从某种模型,在估计 \(I = {EY}\) 时应当尽量利用模型信息。

{\color{red}\textbf{[核心洞察]:} 条件期望法的思想是``用确定性代替随机性''。如果$Y=f(X)+$噪声,直接模拟$Y$要承受噪声的波动;但若能计算$E(Y|X)$,就相当于把噪声``平均掉''了,只留下信号$f(X)$的变化。}

设变量 \(X\) 与 \(Y\) 不独立,根据 Rao-Blackwell 不等式:

\[
\operatorname{Var}\{ \mathrm{E}\left( {Y \mid  X}\right) \}  \leq  \operatorname{Var}\left( Y\right)
\]

{\color{red}\textbf{[为什么方差更小]:} 这来自方差分解公式:$\text{Var}(Y)=E[\text{Var}(Y|X)]+\text{Var}[E(Y|X)]$。条件期望$E(Y|X)$消除了``给定$X$后$Y$的随机波动''$\text{Var}(Y|X)$,只保留了``$X$变化引起的$Y$期望的波动'',因此方差更小。}

又

\[
\mathrm{E}\{ \mathrm{E}\left( {Y \mid  X}\right) \}  = \mathrm{E}Y = I
\]

所以,对 \(Z = \mathrm{E}\left( {Y \mid  X}\right)\) 抽样,用 \(Z\) 的样本平均值来估计 \(I = \mathrm{E}Y\) 比直接用 \(Y\) 的样本平均值的精度更高。这种改善随机模拟估计精度的方法叫做条件期望法, 或 Rao-Blackwell 方法。

例 3.3.5. 设 \(X \sim  {p}_{1}\left( x\right) ,\varepsilon  \sim  \mathrm{N}\left( {0,{\sigma }^{2}}\right)\) 且与 \(X\) 独立,

\[
Y = \psi \left( X\right)  + \varepsilon ,
\]

要估计 \(I = \mathrm{E}Y\) ,可以用条件分布抽样法对二元随机向量 \(\mathbf{Z} = \left( {X,Y}\right)\) 抽样产生 \(Y\) 的样本。从 \(p\left( \cdot \right)\) 抽样得 \({X}_{1},{X}_{2},\ldots ,{X}_{N}\) ,独立地从 \(\mathrm{N}\left( {0,{\sigma }^{2}}\right)\) 抽样得 \({\varepsilon }_{1},{\varepsilon }_{2},\ldots ,{\varepsilon }_{N}\) ,令

\[
{Y}_{i} = \psi \left( {X}_{i}\right)  + {\varepsilon }_{i},i = 1,2,\ldots ,N
\]

然后用 \({Y}_{1},{Y}_{2},\ldots ,{Y}_{N}\) 的样本平均值估计 \({EY}\) :

\[
{\widehat{I}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{Y}_{i}
\]

估计方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  = \frac{\operatorname{Var}\left( {Y}_{1}\right) }{N} = \frac{\operatorname{Var}\left( {\psi \left( {X}_{1}\right) }\right) }{N} + \frac{{\sigma }^{2}}{N}.
\]

另一方面,注意 \(E\left( {Y \mid  X}\right)  = \psi \left( X\right)\) ,也可以只对 \(X\) 抽样然后用条件期望法估计 \({EY}\) :

\[
{\widehat{I}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{i}\psi \left( {X}_{i}\right) ,
\]

估计方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{\operatorname{Var}\left( {\psi \left( {X}_{1}\right) }\right) }{N} < \operatorname{Var}\left( {\widehat{I}}_{1}\right) .
\]

这个例子演示了条件期望法可以缩减误差方差的原因: 对 \(Y\) 抽样分成两步进行: 第一步对 \(X\) 抽样,第二步利用 \(Y \mid  X\) 的条件分布对 \(Y\) 抽样。所以使用 \(Y\) 的样本估计 \({EY}\) 包含了第二步对 \(Y\) 抽样的随机误差,而使用 \(X\) 的函数 \(E\left( {Y \mid  X}\right)  = \psi \left( X\right)\) 的抽样来估计 \({EY}\) 则避免了第二步对 \(Y\) 抽样引起的随机误差。

{\color{red}\textbf{[解释]:} 直接模拟$Y$相当于先抽$X$再加噪声$\varepsilon$——两步都有随机性。条件期望法只抽$X$,用解析公式$\psi(X)$代替加噪声这一步——减少了一个随机源,方差自然变小。这就像用GPS定位代替``GPS+手抖''。}

例 3.3.6. 继续考虑例3.1.1中用随机模拟方法估计 \(\pi\) 的问题。设(X, Y)服从正方形 \(D =\)  \({\left\lbrack  -1,1\right\rbrack  }^{2}\) 上的均匀分布,令 \(\eta  = 1\) 表示(X, Y)落入单位圆 \(C = \left\{  {\left( {x,y}\right)  : {x}^{2} + {y}^{2} \leq  1}\right\}  ,\eta  = 0\) 表示未落入单位圆,则 \(\eta  \sim  \mathrm{B}\left( {1,\pi /4}\right)\) 。设 \(\left( {{X}_{i},{Y}_{i}}\right) ,i = 1,2,\ldots ,N\) 是(X, Y)的独立重复抽样, \({\eta }_{i}\) 表示 \({X}_{i}^{2} + {Y}_{i}^{2} \leq  1\) 是否发生,则 \(\pi\) 的估计为

\[
{\widehat{\pi }}_{1} = \frac{4}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\eta }_{i}
\]

方差为 \(\operatorname{Var}\left( {\widehat{\pi }}_{1}\right)  = \pi \left( {4 - \pi }\right) /N \approx  \frac{2.6968}{N}\) 。

用 \(\zeta  = E\left( {\eta  \mid  X}\right)\) 的样本代替 \(\eta\) 来估计 \({E\eta } = \pi /4\) ,则由

\[
E\left( {\eta  \mid  X = x}\right)  = P\left( {{X}^{2} + {Y}^{2} \leq  1 \mid  X = x}\right)  = P\left( {{Y}^{2} \leq  1 - {x}^{2}}\right)
\]

\[
= P\left( {-\sqrt{1 - {x}^{2}} \leq  Y \leq  \sqrt{1 - {x}^{2}}}\right)
\]

\[
= \sqrt{1 - {x}^{2}}\;\text{ (注意 }Y \sim  \mathrm{U}\left( {-1,1}\right) \text{ ) }
\]

可知 \(\zeta  = \sqrt{1 - {X}^{2}}\) 。其方差为

\[
\operatorname{Var}\left( \zeta \right)  = E\left( {1 - {X}^{2}}\right)  - {\left( E\zeta \right) }^{2} = \frac{2}{3} - \frac{{\pi }^{2}}{16},
\]

估计 \(\pi\) 为

\[
{\widehat{\pi }}_{2} = \frac{4}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\sqrt{1 - {X}_{i}^{2}},
\]

方差为 \(\operatorname{Var}\left( {\widehat{\pi }}_{2}\right)  \approx  {0.7971}/N,\operatorname{Var}\left( {\widehat{\pi }}_{1}\right) /\operatorname{Var}\left( {\widehat{\pi }}_{2}\right)  \approx  {3.4}\) 。

另外, \(\zeta\) 仅依赖于 \({X}^{2}\) ,容易发现若 \(U \sim  \mathrm{U}\left( {0,1}\right)\) 则 \({X}^{2}\) 和 \({U}^{2}\) 同分布,所以可取 \(\xi  = \sqrt{1 - {U}^{2}}\) ,其中 \(U \sim  \mathrm{U}\left( {0,1}\right)\) 。这时,函数 \(h\left( u\right)  = \sqrt{1 - {u}^{2}}\) 是 \(u \in  \left( {0,1}\right)\) 的单调函数,可以利用对立变量法,构造 \(\pi\) 的估计量为

\[
{\widehat{\pi }}_{3} = \frac{4}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{\sqrt{1 - {U}_{i}^{2}} + \sqrt{1 - {\left( 1 - {U}_{i}\right) }^{2}}}{2},
\]

其中 \({U}_{1},{U}_{2},\ldots ,{U}_{N}\) 为 \(U\left( {0,1}\right)\) 随机数,则 \(\operatorname{Var}\left( {\widehat{\pi }}_{3}\right)  \approx  {0.11}/N,\operatorname{Var}\left( {\widehat{\pi }}_{1}\right) /\operatorname{Var}\left( {\widehat{\pi }}_{3}\right)  \approx  {25}\) 。

也可以用对立变量法来改进 \({\widehat{\pi }}_{2}\) 。令 \(W = {U}^{2} - \frac{1}{3}\) ,则 \({EW} = 0\) 且 \(W\) 与 \(\zeta\) 负相关。可以先进行一个小规模的模拟估计 \(\operatorname{Cov}\left( {\zeta ,W}\right)\) 和 \(\operatorname{Var}\left( W\right)\) 得到 \(b =  - \operatorname{Cov}\left( {\zeta ,W}\right) /\operatorname{Var}\left( W\right)\) 的近似值,用 \(\zeta \left( b\right)  = \zeta  + {bW}\) 代替 \(\zeta\) 进行抽样,可以减小 \({\widehat{\pi }}_{2}\) 的方差。

例 3.3.7. 设随机变量 \(\mathbf{Y} \sim  f\left( \mathbf{y}\right)\) ,为了估计 \(\theta  = {Eh}\left( \mathbf{Y}\right)\) ,经常使用标准化重要抽样法

\[
{\widehat{\theta }}_{1} = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}h\left( {\mathbf{X}}_{i}\right) }{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}}, \tag{3.55}
\]

其中 \({\mathbf{X}}_{i}\) 为试抽样密度 \(g\left( \mathbf{x}\right)\) 的样本,权重 \({W}_{i} = f\left( {\mathbf{X}}_{i}\right) /g\left( {\mathbf{X}}_{i}\right)\) 。如果 \(\mathbf{x} = \left( {\mathbf{u},\mathbf{v}}\right) ,h\left( \mathbf{x}\right)  =\)  \({h}_{1}\left( \mathbf{u}\right)\) ,则只要对 \(\mathbf{X} = \left( {\mathbf{U},\mathbf{V}}\right)\) 的分量 \(\mathbf{U}\) 抽样得到 \({\mathbf{U}}_{i},i = 1,\ldots ,N\) 及新的权重 \({\widetilde{W}}_{i} =\)  \({f}_{\mathbf{U}}\left( {\mathbf{U}}_{i}\right) /{g}_{\mathbf{U}}\left( {\mathbf{U}}_{i}\right)\) ,其中 \({f}_{\mathbf{U}}\left( \mathbf{u}\right)  = \int f\left( {\mathbf{u},\mathbf{v}}\right) d\mathbf{v},{g}_{\mathbf{U}}\left( \mathbf{u}\right)  = \int g\left( {\mathbf{u},\mathbf{v}}\right) d\mathbf{v}\) ,则可估计 \(\theta\) 为

\[
{\widehat{\theta }}_{2} = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{\widetilde{W}}_{i}{h}_{1}\left( {\mathbf{U}}_{i}\right) }{\mathop{\sum }\limits_{{i = 1}}^{N}{\widetilde{W}}_{i}}, \tag{3.56}
\]

这时

\[
\operatorname{Var}\left( {\widetilde{W}}_{i}\right)  \leq  \operatorname{Var}\left( {W}_{i}\right) ,i = 1,2,\ldots ,N. \tag{3.57}
\]

事实上,

\[
\frac{{f}_{\mathbf{U}}\left( \mathbf{u}\right) }{{g}_{\mathbf{U}}\left( \mathbf{u}\right) } = \int \frac{f\left( {\mathbf{u},\mathbf{v}}\right) }{{g}_{\mathbf{U}}\left( \mathbf{u}\right) }{dv} = \int \frac{f\left( {\mathbf{u},\mathbf{v}}\right) }{{g}_{\mathbf{U}}\left( \mathbf{u}\right) {g}_{\mathbf{V} \mid  \mathbf{U}}\left( {\mathbf{v} \mid  \mathbf{u}}\right) }{g}_{\mathbf{V} \mid  \mathbf{U}}\left( {\mathbf{v} \mid  \mathbf{u}}\right) d\mathbf{v}
\]

\[
= \int \frac{f\left( {\mathbf{u},\mathbf{v}}\right) }{g\left( {\mathbf{u},\mathbf{v}}\right) }{g}_{\mathbf{V} \mid  \mathbf{U}}\left( {\mathbf{v} \mid  \mathbf{u}}\right) d\mathbf{v} = {E}_{g}\left\{  {\left. \frac{f\left( {\mathbf{U},\mathbf{V}}\right) }{g\left( {\mathbf{U},\mathbf{V}}\right) }\right| \;\mathbf{U} = \mathbf{u}}\right\}   \tag{3.58}
\]

于是

\[
{\operatorname{Var}}_{g}\left\{  \frac{f\left( {\mathbf{U},\mathbf{V}}\right) }{g\left( {\mathbf{U},\mathbf{V}}\right) }\right\}   \geq  {\operatorname{Var}}_{g}\left\{  {{E}_{g}\left\lbrack  {\left. \frac{f\left( {\mathbf{U},\mathbf{V}}\right) }{g\left( {\mathbf{U},\mathbf{V}}\right) }\right| \;\mathbf{U}}\right\rbrack  }\right\}   = {\operatorname{Var}}_{g}\left\{  \frac{{f}_{\mathbf{U}}\left( \mathbf{U}\right) }{{g}_{\mathbf{U}}\left( \mathbf{U}\right) }\right\}  . \tag{3.59}
\]

\subsection{随机数复用}

在统计研究中, 经常需要比较若干种统计方法的性能, 如偏差、方差、覆盖率等。除了努力获取理论结果以外,可以用随机模拟方法进行比较: 重复生成 \(N\) 组随机样本,对每个样本同时用不同统计方法计算结果,最后从 \(N\) 组结果比较不同方法的性能。这样比较时,并没有对每种方法单独生成 \(N\) 组样本,而是每个样本同时应用所有要比较的方法。这样不仅减少了计算量,而且在比较时具有更高的精度。

{\color{red}\textbf{[为什么复用更精确]:} 这是``配对比较''思想——就像医学试验中同一病人服药前后对比,比两组不同病人对比更可靠。用相同随机数时,两种方法的误差高度相关,相减时大部分误差抵消,差值的方差远小于独立样本的情况。}

例 3.3.8. 对正态分布总体 \(X \sim  N\left( {\mu ,{\sigma }^{2}}\right)\) ,如果有样本 \({X}_{1},{X}_{2},\ldots ,{X}_{n}\) ,估计 \({\sigma }^{2}\) 有两种不同的公式:

\[
{\widehat{\sigma }}_{1}^{2} = \frac{1}{n - 1}\mathop{\sum }\limits_{{i = 1}}^{n}{\left( {X}_{i} - \bar{X}\right) }^{2},{\widehat{\sigma }}_{2}^{2} = \frac{1}{n}\mathop{\sum }\limits_{{i = 1}}^{n}{\left( {X}_{i} - \bar{X}\right) }^{2}. \tag{3.60}
\]

希望比较两个估计量的偏差和均方误差:

\[
{b}_{1} = E{\widehat{\sigma }}_{1}^{2} - {\sigma }^{2},\;{b}_{2} = E{\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}, \tag{3.61}
\]

\[
{s}_{1} = E{\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2},\;{s}_{2} = E{\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}. \tag{3.62}
\]

当然, 这个问题很简单, 可以得到偏差和均方误差的理论值:

\[
{b}_{1} = 0,\;{b}_{2} =  - \frac{1}{n}{\sigma }^{2}, \tag{3.63}
\]

\[
{s}_{1} = \frac{2{\sigma }^{4}}{n - 1},\;{s}_{2} = \frac{\left( {{2n} - 1}\right) {\sigma }^{4}}{{n}^{2}},\;{s}_{1} - {s}_{2} = \frac{\left( {{3n} - 1}\right) {\sigma }^{4}}{{n}^{2}\left( {n - 1}\right) }. \tag{3.64}
\]

我们用随机模拟来作比较。重复地生成 \(N\) 组样本 \(\left( {{X}_{1}^{\left( j\right) },{X}_{2}^{\left( j\right) },\ldots ,{X}_{n}^{\left( j\right) }}\right) ,j = 1,2,\ldots ,N\) 。对每组样本分别计算 \({\widehat{\sigma }}_{1j}^{2} = \frac{1}{n - 1}\mathop{\sum }\limits_{{i = 1}}^{n}{\left( {X}_{i}^{\left( j\right) } - {\bar{X}}^{\left( j\right) }\right) }^{2}\) 和 \({\widehat{\sigma }}_{2j}^{2} = \frac{1}{n}\mathop{\sum }\limits_{{i = 1}}^{n}{\left( {X}_{i}^{\left( j\right) } - {\bar{X}}^{\left( j\right) }\right) }^{2}\) ,得到偏差和均方误差的估计

\[
{\widehat{b}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{\widehat{\sigma }}_{1j}^{2} - {\sigma }^{2},\;{\widehat{b}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{\widehat{\sigma }}_{2j}^{2} - {\sigma }^{2}, \tag{3.65}
\]

\[
{\widehat{s}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{\left( {\widehat{\sigma }}_{1j}^{2} - {\sigma }^{2}\right) }^{2},\;{\widehat{s}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{\left( {\widehat{\sigma }}_{2j}^{2} - {\sigma }^{2}\right) }^{2}. \tag{3.66}
\]

容易看出, 两种方法使用相同的模拟样本得到的偏差、均方误差的估计精度与每种方法单独生成模拟样本得到的估计精度相同。

但是,如果要估计 \({\Delta s} = {s}_{1} - {s}_{2}\) ,利用相同的样本的估计精度更好。一般地,

\[
\operatorname{Var}\left( {{\widehat{s}}_{1} - {\widehat{s}}_{2}}\right)  = \operatorname{Var}\left( {\widehat{s}}_{1}\right)  + \operatorname{Var}\left( {\widehat{s}}_{2}\right)  - 2\operatorname{Cov}\left( {{\widehat{s}}_{1},{\widehat{s}}_{2}}\right) . \tag{3.67}
\]

{\color{red}\textbf{[关键项]:} 注意这个$-2\text{Cov}(\hat{s}_1,\hat{s}_2)$项!当两种方法用相同样本时,协方差通常很大且为正(因为都依赖同一数据),这一项能大幅降低方差。这就是``共同误差源相互抵消''的数学体现。}

如果每种方法使用不同的样本, 则 (3.67) 变成

\[
\operatorname{Var}\left( {{\widehat{s}}_{1} - {\widehat{s}}_{2}}\right)  = \frac{1}{N}\left\lbrack  {\operatorname{Var}\left( {\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2}\right)  + \operatorname{Var}\left( {\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}\right) }\right\rbrack  . \tag{3.68}
\]

如果两种方法利用相同的样本计算, 则 (3.67) 变成

\[
\operatorname{Var}\left( {{\widehat{s}}_{1} - {\widehat{s}}_{2}}\right)  = \frac{1}{N}\left\lbrack  {\operatorname{Var}\left( {\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2}\right)  + \operatorname{Var}\left( {\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}\right)  - 2\operatorname{Cov}\left( {{\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2},{\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}}\right) }\right\rbrack  .
\]

(3.69)

而 \({\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2}\) 与 \({\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}\) 明显具有强正相关 (见习题10),所以两种方法针对相同样本计算时 \({\widehat{s}}_{1} - {\widehat{s}}_{2}\) 的方差比两种方法使用单独样本时的方差要小得多。这说明重复利用相同的随机数或样本往往可以提高比较的精度。

{\color{red}\textbf{[实用建议]:} 随机数复用是最容易实施的方差缩减技巧——几乎不需要数学分析,只需在比较算法时``用同一组随机数''即可。这应该成为所有模拟比较研究的标准做法。}
